DOE/ER/25053-5 
UC-32 


Courant  Institute  of 
Mathematical  Sciences 


Laplace's  Equation  and  the  Dirichlet-Neumann  Map 
in  Multiply  Connected  Domains 

Anne  Greenbaum,  Leslie  Greengard  and 
Geoffrey  B.  McFadden 


iupported  by  U.S.  Department  of  Energy 
jrant  No.  DE-FGO2-88ER25053 


m     -H  S  o  vlathematics  and  Computers 

O    (l»   4_l    g    d)  ' 

uT)  c  CO  a  c  Vlarch  1991 

tN  c  3  a;  c 
^  <  CPS  o 

;\  s  w  <u  >, 
o  (a  <D  jz  cu 

C    (0  -H  JJ 

0)    nH     ^   rH 

>H    1-1    (C   Q    E 

S  O  i-l 


NEW  YORK  UNIVERSITY 


UNCLASSIFIED 

DOE/ER/25053-5 

UC-32 

Mathematics  and  Computers 


Courant  Institute  of  Mathematical  Sciences 
New  York  University 


Laplace's  Equation  and  the  Dirichlet-Neumann  Map 
in  Multiply  Connected  Domains 

A.  Greenbaum,  L.  Greengard,  and  G.  B.  McFadden 
March  1991 


The  first  two  authors  were  supported  by  U.  S.  Department  of  Energy  Grant  DE-FG02- 
8SER25053.  The  third  author  is  at  the  Compiiting  and  Applied  Mathematics  Laboratory, 
National  Institute  of  Standards  and  Technology.  He  was  supported  by  the  NASA  Micro- 
gravity  Science  and  Applications  Program,  and  the  DARPA  Applied  and  Computational 
Mathematics  Program. 

UNCLASSIFIED 


DISCLAEVIER 


This  report  was  prepared  as  an  account  of  work  sponsored  by  an  agency  of  the 
United  States  Government.  Neither  the  United  States  Government  nor  any  agency 
thereof,  nor  any  of  their  employees,  makes  any  warranty,  express  or  implied,  or  as- 
sumes any  legal  liability  or  responsibihty  for  the  accuracy,  completeness,  or  useful- 
ness of  any  information,  apparatus,  product,  or  process  disclosed,  or  represents  that 
its  use  would  not  infringe  privately  owned  rights.  Reference  herein  to  any  specific 
commerciaJ  product,  process,  or  service  by  trade  name,  trademark,  manufacturer, 
or  otherwise,  does  not  necessarily  constitute  or  imply  its  endorsement,  recommen- 
dation, or  favoring  by  the  United  States  Government  or  any  agency  thereof.  The 
views  and  opinions  of  authors  expressed  herein  do  not  necessarily  state  or  reflect 
those  of  the  United  States  Government  or  any  agency  thereof. 


Printed  in  U.S.A. 

Available  from 

National  Technical  Information  Service 

U.S.  Department  of  Commerce 

5285  Port  Royal  Road 

Springfield,  VA  22161 


11 


Table  of  Contents 

Page 
Abstract  1 

1.  Introduction  2 

2.  Siniply  Connected  Domains  3 

3.  Multiply  Connected  Domains  5 

3.1  A  direct  formulation  6 

3.2  The  exterior  problem      '  7 

4.  Formulation  of  Numerical  Method  8 

4.1  Solution  of  the  Discrete  System  10 

5.  The  Dirichlet-Neumann  Map  and  the  Neumann  Problem  12 

5.1  The  Neumann  Problem  13 

6.  Numerical  Results  14 

7.  Conclusions  17 
References  18 
Figures  21 


111 


Laplace's  Equation  and  the 

Dirichlet-Neumann  Map  in  Multiply 

Connected  Domains 

Anne  Greenbaum*  Leslie  Greengard^ 

Geoffrey  B.  McFadden* 

March,  1991 

Abstract 

A  variety  of  problems  in  material  science  and  fluid  dynamics  require 
the  solution  of  Laplace's  equation  in  multiply  connected  domains.  In- 
tegral equation  methods  are  natural  candidates  for  such  problems,  since 
they  discretize  the  boundary  alone,  require  no  special  effort  for  free  bound- 
aries, and  achieve  superalgebraic  convergence  rates  on  sufficiently  smooth 
domains  in  two  space  dimensions,  regardless  of  shape.  Current  integral 
equation  methods  for  the  Dirichlet  problem,  however,  require  the  solution 
of  A/  independent  problems  of  dimension  A',  where  M  is  the  number  of 
boundary  components  and  A'  is  the  total  number  of  points  in  the  dis- 
cretization. In  this  paper,  we  present  a  new  boundary  integral  equation 
approach,  valid  for  both  interior  and  exterior  problems,  which  requires 
the  solution  of  a  single  linear  system  of  dimension  A'  +  M.  We  solve  this 
system  by  making  use  of  an  iterative  method  (GMRES)  combined  willi 
the  fast  multipole  method  for  the  rapid  calculation  of  the  necessary  matrix 
vector  products.  For  a  two-dimensional  system  with  200  components  and 
100  points  on  each  boundary,  we  gain  a  speedup  of  a  factor  of  100  from 
the  new  analytic  formulation  and  a  factor  of  50  from  the  fast  multipole 
method.  The  resulting  scheme  brings  large  scale  calculations  in  extremely 
complex  domains  within  practical  reach. 
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2  A.  Greenbaum,  L.  Greengard  and  G.  B.  McFadden 

1      Introduction 

Several  problems  in  applied  mathematics  involve  the  solution  of  Laplace's  equa- 
tion in  multiply  connected  domains.  Examples  arise  in  fluid  dynamics,  electro- 
statics, and  various  aspects  of  material  science  (see,  for  example,  [2],  [6],  [12], 
[14],  [16],  [18],  [19],  [26],  [27]).  In  some  circumstances,  the  solution  is  desired  at 
a  set  of  points  in  the  interior  of  the  domain.  In  other  circumstances,  values  of 
the  normal  derivative  are  to  be  derived  from  a  specification  of  function  values  on 
the  boundary  (the  Dirichlet-Neumann  map).  In  the  numerical  treatment  of  such 
problems,  second  kind  integral  equation  methods  are  a  natural  choice,  having  a 
number  of  advantages  over  finite  difference  or  finite  element  discretizations.  The 
latter  require  a  discretization  of  the  entire  computational  domain,  arc  difficult 
to  apply  to  exterior  problems,  lead  to  poorly  conditioned  linear  systems,  and 
achieve  a  low  order  of  accuracy,  although  the  resulting  matrices  are  sparse  and 
amenable  to  efficient  solution  by,  for  example,  a  multigrid  procedure.  Second 
kind  integral  equation  methods  require  a  discretization  of  the  boundary  alone 
(so  that  the  size  of  the  resulting  linear  system  is  drastically  reduced),  achieve 
superalgebraic  convergence  on  smooth  two-dimensional  domains,  are  as  easily 
applied  to  exterior  as  to  interior  problems,  and  are  well  condil  ioned.  They  do, 
however,  give  rise  to  dense  matrices. 

In  simply  connected  domains,  the  most  obvious  method  for  solving  the  linear 
systems  which  result  from  discretizing  integral  equations  is  (Jaussian  elimina- 
tion, which  requires  O(N^)  operations  where  A^  is  the  number  of  points  on  the 
boundary.  On  the  other  hand,  a  number  of  authors  have  observed  ihat  since  such 
systems  are  well-conditioned,  iterative  methods  are  easily  applied  and  achieve 
a  fixed  accuracy  in  a  number  of  iterations  independent  of  A'  (sec,  for  example, 
[2],  [7],  [21]).  Each  iteration,  of  course,  would  appear  to  re(|uire  0{N^)  opera- 
tions so  that  the  net  cost  of  solving  the  linear  system  would  he  0{N').  During 
the  past  few  years,  however,  several  schemes  have  appeared  which  reduce  the 
cost  of  each  iteration  to  O(A'logA')  or  0{N)  [1,4,5,10,11,20,21,25],  and  their 
application  results  in  optimal  order  methods  for  solving  Laplace's  equation  Of 
particular  note  is  the  paper  of  Rokhlin  [21],  in  which  such  a  hybrid  method  first 
appeared. 

In  multiply  connected  domains,  unfortunately,  the  sitiialion  is  not  so  sim- 
ple. Second  kind  integral  equation  methods  for  the  Diriclilct  i^roblcm  develop  a 
nullspace  of  dimension  equal  to  the  number  of  boundary  components  A/ .  In  or- 
der to  obtain  a  solution,  therefore,  one  of  two  strategies  has  been  employed.  The 
first  is  to  project  the  Dirichlet  data  onto  the  range  of  the  operator  by  subjecting 
the  data  to  certain  compatibility  constraints  (see,  for  example,  [3],  [IC]).  While 
the  algorithms  of  [3]  and  [16]  differ  in  several  details,  both  re(|uire  the  solution 
of  A/  adjoml  integral  equations,  each  of  dimension  A',  where  .^  is  the  number 
of  points  in  the  discretization  of  the  entire  boundary.  An  alienialive  approach, 
due  to  Mayo  [15],  uses  a  modified  integral  equation  descrilu-d  i)y  Mikhliii  [18], 
This  scheme  does  not  require  solutions  of  the  adjoint  equation,  but  does  require 
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the  solution  of  M  (non-singular)  integral  equations  of  dimension  A'.  With  the 
use  of  the  fast  summation  techniques  mentioned  above,  the  work  reciuired  by 
the  above  methods  is  of  the  order  0(N  ■  A/). 

In  this  paper,  we  present  an  algorithm  for  solving  Laplace's  equation  in  M- 
ply  connected  domains  in  two  space  dimensions  where  M  is  large  (up  to  several 
hundred),  in  which  case  the  above  schemes  undergo  a  serious  degradation  in 
performance.  The  asymptotic  CPU  time  requirements  of  our  algorithm  are  of 
the  order  0{N  +  M),  0{N  +  M')  or  0(N  +  AP),  depending  on  the  details 
of  the  implementation,  rather  than  0{N  ■  M).  For  a  system  with  200  compo- 
nents and  100  points  on  each  boundary  (M  =  200,  A^  =  20,000),  we  gain  a 
speedup  of  a  factor  of  100  from  the  new  analytic  formulation.  Another  factor 
of  50  in  speedup  is  gained  from  using  the  fast  multipole  method  [5,10,21]  to 
apply  the  discretized  integral  operators.  The  resulting  scheme  brings  large  scale 
calculations  in  extremely  complex  domains  within  practical  reach. 

The  paper  is  organized  cis  follows:  Section  2  reviews  the  integral  equation 
approach  to  the  Dirichlet  problem  in  simply  connected  domains.  Section  3  ex- 
tends the  analysis  to  multiply  connected  domains,  and  Section  4  describes  our 
numerical  method  in  some  detail.  Section  5  explains  how  to  obtain  the  Dirichlet- 
Neumann  map  and  how  to  solve  the  Neumann  problem.  Section  6  pre.seiits  the 
results  of  some  numerical  experiments,  and  Section  7  contains  our  conclusions. 

2      Simply  Connected  Domains 

We  will  briefly  describe  the  solution  of  the  Dirichlet  problem  for  Laplace's  equa- 
tion in  simply  connected  domains  via  classical  potential  theory.  Descript  ions  of 
this  theory  may  be  found,  for  example,  in  [8],  [9],  [13],  and  [18]. 

Let  us  begin  by  considering  a  finite  open  region  D  in  the  plane  willi  liouiidary 
L,  which  we  cissume  to  be  smooth  and  to  have  continuous  curvature.  We  denote 
by  E  the  open  region  in  the  plane  exterior  to  L.  The  interior  Dirichlet  problem, 
then,  involves  the  determination  of  a  function  U(P)  which  satisfies 

A(/(F)  =  0         ioT  P  e  D  (1) 

and  the  boundary  condition 

I'm  U{P)  =  f{Q)         for  Q€  L  .  (2) 

P  —  Q 
P€0 


We  seek  the  solution  U{P)  in  the  form  of  a  double  layer  poteiitia 

'  di'Q 

where  P  is  a  point  inside  D.  ^'{Q)  is  the  value  of  the  unknown  dipole  dist  ril)ution 
at  the  boundary  point  Q.  and  ^^  represents  the  outward  normal  derivative  at 
the  point  Q 


U[P)^^    f  n(Q)i-\u\Q-P\dQ  ,  (3) 
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For  any  point  Qo  on  the  boundary,  the  double  layer  potential  satisfies  the 
jump  relations 


^/i(Qo)  +^  f  t^(Q)j^  In  \Q  -  Qo\  dQ  (4) 

2  2n  J I  duQ 


lim    U{P) 

P—Qo 

PSD 

lim    f/(P)     =     -LiQ,)+^   f  ^,(Q)^\n\Q-Qo\dQ.  (5) 

^^^0  /  Z7r   If 

PiE  •' '' 


duQ 


In  order  to  satisfy  the  boundary  condition  (2),  the  relation  (-1)  implies  that  the 
function  /i  must  satisfy  the  integral  equation 


u. 


^(Qo)+-  /  A'(Q)^— ln|g-QoMg  =  2/(Qo) 


L  di'Q 


(6) 


We  note  that  the  kernel  g^  In  |Q  —  Qo\  is  continuous  along  the  curve;  in  fact, 
for  future  reference, 


lim    J-\n\Q-Qo\=\K(Qo], 

at:  L  V 


(7) 


where  k(Q)  denotes  the  curvature  of  L  at  the  point  Q.  Moreover,  it  is  well  known 
that  equation  (6)  has  no  nontrivial  homogeneous  solutions.  W'c  may,  therefore, 
conclude  from  the  Fredholm  alternative  that  (6)  has  a  uni(|U(>  solution  for  any 
integrable  data  f{Q). 

For  the  exterior  Dirichlet  problem. 


AU{P)     =     0         (or  P  e  E 
lirn  U(P)     =     !(Q)         for  Q  €  Z,  , 


P  —  Q 


(8) 


we  encounter  an  obvious  difficulty  if  we  seek  the  solution  in  terms  of  a  double 
layer  potential  alone.  Such  a  representation  decays  to  zero  al  infinity,  while  an 
appropriate  condition  for  uniqueness  is  that  U  tend  to  a  constant  If  we  were 
to  try  to  use  the  same  double  layer  potential  as  before,  the  jump  relation  (5) 
would  yield  the  integral  equation 


-^^(Q 


o)+-   f  fi(Q)^\n\Q-Qo\dQ='2f(Q,) 


(9) 


The  Fredholm  alternative  still  applies,  of  course,  but  the  homogeneous  equation 
has  the  nontrivial  solution  ^i  =  C,  &  constant.  Let  us,  therefore,  seek  a  slightly 
different  representation,  namely 

r(P)     =     ^   hi(Q)^\n\Q-P\dQ+-^  I  ,,(Q)dQ 


Q) 


duq 


ln|Q-P|+  1 


dQ 


;io) 
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The  integral  equation  resulting  from  application  of  the  jump  relation  (5)  is  then 

d 


u. 


^l{Qo)  +  -  /  ^i{Q) 


^,^^in|g-go|  +  i 


dQ  =  2f(Qo).  (II) 


A  proof  of  the  fact  that  (11)  has  a  unique  solution  for  any  integrable  data  f{Q) 
can  be  found  in  [18]. 

3      Multiply  Connected  Domains 

Let  us  now  consider  the  Dirichlet  problem  in  a  finite  open  region  D  in  the  plane 
which  is  (M  +  l)-ply  connected  (Fig.  1).  The  outer  boundary  of  D  will  be  de- 
noted by  Lo,  while  the  interior  boundary  curves  will  be  denoted  by  Li ,  .  .. ,  La/  . 
If  we  seek  the  solution  U(P)  as  a  double  layer  potential,  and  u.se  llic  jump 
relation  (4),  we  obtain  the  integral  equation 

fi{Qo)+-  /  f^{Q)-r—\n\Q  -  Qo\dQ  =  2f(Qo)  ,  (12) 

where  f(Q)  is  the  given  boundary  data. 

Remark  3.1:  We  will  follow  the  standard  convention  that  for  i;i/f  no;  problems, 
Q—  refers  to  the  normal  derivative  in  the  direction  outward  from  the  domain 
D.  Thus,  whether  the  boundary  point  Qq  lies  on  the  outer  boundary  or  on  one 
of  the  interior  curves,  the  relevant  jump  condition  is  (4).  For  er/f  nor  problems, 
g—  refers  to  the  inward  normal  derivative  and  the  relevant  jump  condition  is 
(5). 

There  is  a  fundamental  problem  with  the  linear  system  (12),  wliirli  is  that 
there  are  M  nonlrivial  homogeneous  solutions  of  the  form 

r    1     forQGL.  . 

That  they  are  null  solutions  follows  immediately  from  the  identity 


hi 


c,  r    1      for  Qo  inside  Z-, 

■^\n\Q-Qo\dQ=  }    i     for  Qo  on  L, 

^  [    0      for  Qo  outside  L,. 


For  a  proof  that  they  span  the  nullspace,  see  [8]. 

In  short,  wc  may  conclude  from  Fredholm  theory  that  equation  (12)  can  be 
solved  only  if  the  right-hand  side  is  orthogonal  to  the  A/  independent  solutions 
of  the  adjomt  mlegral  equation 

This  approach  is  the  one  taken  in  [3]  and  [16],  and  requires  the  dci(-rniiiiation 
of  the  A/  independent  solutions  of  the  adjoint  system. 


<T(Qo)+  -  /t(Q)^^^  In  |Q-Qo|(fQ=0 
^  Jl 
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3.1      A  direct  formulation 

A  somewhat  different  approach  to  the  multiply-connected  probleni  is  suggested 
by  Mikhlin  [18].  Rather  than  making  use  of  the  general  Fredhohn  theory,  as  we 
did  above,  let  us  seek  a  solution  in  the  form 


d 


M 


(13) 


i  =  l 


where  St  is  a  point  inside  the  interior  boundary  curve  Lt-.  In  the  context  of 
fluid  dynamics,  the  unknown  constants  At  correspond  to  circulation,  but  we 
will  think  of  them  as  playing  the  role  of  Lagrange  multipliers.  Applying  the 
jump  relations  as  before,  the  resulting  integral  equation  is 


^(Qo)  +  -  / 

2/(Qo)-2^/4iln|Qo-5, 


^i{Q)-^\u\Q-Qo\dQ  = 


M 


14) 


i  =  l 


As  an  integral  equation  for  the  dipole  distribution  /u,  the  rank  deficiency 
encountered  previously  is  still  present.  We  have  only  modified  the  right-hand 
side.  The  next  step,  however,  is  to  replace  equation  (14)  by  the  following; 


^  Jl 


iQ) 


^— ln|Q-Qo|  +  5(Q,Qo) 

OUq 

M 


dQ  = 


2/(Qo)-2^/l,ln|P-.9,| 


(15) 


k  =  i 


whe 


1     ifQ,Qo6Lt,<-=  1, 
0     otherwise. 


A/ 


S(Q,Qo)  = 
This  modification  of  the  kernel  adds  a  term  of  the  form 

-  /  ^'{Q)dQ 

to  equation  (14)  if  and  only  if  Qo  lies  on  an  interior  boundary  curve  Lj.  The 
value  of  this  modification  is  that  the  integral  equation  (15)  is  nonsingular  -  it 
has  a  unique  solution  for  any  integrable  right-hand  side.  A  proof  may  be  found 
in  [18]. 

Notice,  however,  that  we  have  made  no  use  as  yet  of  the  unknown  constants 
Ak  Given  these  A/  extra  degrees  of  freedom,  we  may  subject  the  system  (15) 
to  the  constraint  conditions 


L 


,i(Q)dQ  =  0 


k=  1 A/ 


(16) 
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But  then  (15)  and  (14)  are  identical,  so  that  we  have  constructed  a  solution  to 
our  problem  in  the  desired  form  (13). 

Let  us,  therefore,  consider  a  single  system  of  equations,  namely  (M)  and 
(16)  taken  together  and  rewritten  in  the  form 


/^(Qo)+  -  /  fi(Q)-^—\n\Q-Qo\dQ  + 

T!  J I  dUQ 

M 

2^/ltln|go-5i|  =  2/(Qo) 


i  =  l 


L 


(17) 


^i(Q)dQ  =  Q 


1,...,A/  . 


By  the  preceding  discussion,  this  system  is  invertible  and  simultaneously  deter- 
mines the  values  of  the  constant  At  and  the  desired  dipole  density  ii(Q). 

3.2      The  exterior  problem 


Let  E  be  the  open  region  in  the  plane  exterior  to  M  closed  curves  A i, 
(Fig    2).  As  suggested  by  Mikhlin  [18],  we  seek  a  solution  in  the  form 

f^(P)=i-   f  ^(Q)^\u\Q-P\dQ+l-   f  M{Q)dQ 
2t  Jl  0//Q  27r  y^ 


'A/ 


M 


-^-Y^AtlulP-Stl  , 


where  St  is  a  point  inside  the  boundary  curve  Lt-  We  will  require  thai 

M 

Y,A,=0 
t  =  i 

to  ensure  that  ('  is  bounded  at  infinity. 

Application  of  the  jump  relation  (5)  yields  the  integral  equation 


(18) 


(19) 


;/ 


-^i(Qo)+-   I   fj(Q)-^\n\Q-Qo\dQ+-   j   ^(Q)dQ^ 


M 

2J{Qa)-2Y^A,\n\Qo-S,\ 


(20) 


k^\ 


whicli  ha-s  M  —  1  nontriviaj  homogeneous  solutions.  This  equation  can  br  mod- 
ified ci-s  above,  rcplacmg  equation  (20)  with 


-/'(C^n)  + 


i/. 


^i{Q) 


d 
duQ 


ln|Q-Qo|  +  /?(Q,Qo) 


dQ  + 
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1    [  " 

-  /  /i(Q)rfQ  =  2/(Qo)-2^/l,ln|Qo-5,|, 


(21 


whe 


R(Q,Qo)=\\    i<;^^^oGL,,^=l M 

0     otherwise. 


Note  that  the  A/""  boundary  curve  is  distinguished  from  the  others  in  the  def- 
inition of  R{Q,Qo)-  The  value  of  this  modification  of  the  kernel  is  that,  as  for 
the  interior  problem,  the  resulting  system  (21)  is  invertible  (see  [18]).  We  now 
make  use  of  the  unknown  constants  At-  Since  we  have  already  required  that 
they  satisfy  (19),  we  are  left  with  A/  -  1  degrees  of  freedom.  \VV  may,  therefore, 
subject  the  preceding  linear  system  to  the  constraints 


/. 


fi{Q)dQ  =  0        k  =  l M  -\  .  (22) 

Equations  (20)  and  (21)  are  then  identical,  and  we  have  constructed  a  solution 
to  the  exterior  problem  in  the  form  (18).  Rewriting  (19),  (20),  and  (22)  as  a 
simultaneous  set  of  equations  for  the  unknown  constants  Ai,  and  tiie  unknown 
dipole  density  /i,  we  have 


u. 


f^iQo)--  /  /i(Q) 


^ln|Q-goI+l 


dQ 


M 


-2^>ltln|Qo-5t|  =  -2/(Qo) 


M 


I  /i(Q) 


Y.At  =  Q  (23) 

k  =  \ 

dQ  =  0        k=  I A/-  1  . 


4      Formulation  of  Numerical  Method 

In  order  to  solve  the  systems  (17)  and  (23),  we  use  the  Nyslrbiii  algorithm  ba.scd 
on  the  trapezoidal  rule.  We  use  the  trapezoidal  rule  for  (luadralure  since  it 
achieves  superalgebraic  convergence  for  smooth  function.s  on  periodic  domains. 
In  more  detail,  we  assume  that  we  are  given  Nt  points  equispaced  in  arciength 
on  each  boundary  component  Lt  and  associate  with  each  siicli  point,  denoted 
Q,  ,  an  unknown  density  value  /jf .  The  step  in  arciength  in  the  discretization 
of  Lk  will  be  denoted  /it  =  \Lt\/Nt,  where  \Lt\  is  the  length  of  the  curve  Lt. 
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The  total  number  of  discretization  points  is  denoted  by 


N     -- 

--  E^* 

for  the  exterior  problem 

N     -- 

M 

for  the  interior  problem 

ifc  =  0 


We  then  replace  the  equations  in  (17)  by 


/^:  +  iE''*i:^'ni^;^-^:i/^.^  + 


k  =  0  7  =  1  V, 

A/ 


2^^,ln|g;-5,|   =  2f(Q[) 


t=i 


Nk 


and  the  equations  in  (23)  by 


M  N» 


'■.-iZX'Z 


k=]  ;=1 


duQ. 


iniQ;-Q;i  +  i 


^'.  - 


M 


2Y^A,\n\Q\-St\   =    -2f{Q[] 


(24) 


M 

t=i 


Ak  =  Q 


(25) 


^/j*/ii  =0         k  =  l M-1, 

j  =  i 

Remark  4.1:   In  the  preceding  linear  systems,  the  diagonal  terms  of  the  form 
3^1n|Q^*^  -  Q',\  should  be  replaced  by  the  appropriate  limit  k(Q\)/'1  when 

1*  '_  ni 


Q    =  Q[  in  accordance  with  equation  (7). 
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The  discrete  equations  for  the  interior  and  exterior  problems  may  be  written 
respectively  as 


C         D-   J    [a  J    -    [    0 


(26) 
(27) 


where,  for  the  interior  problem, 


P=(f^° /i°N„,/i},.../i)v, /if',-     ,/'Jlf„) 

is  the  vector  of  unknown  density  values, 

is  the  vector  of  unknown  coefficients,  and 

/  =  (/] .  •  •  • .  Ino '  /i .  ■  ■  ■ .  /n,  '  •  ■  ■ .  /i    '  •  •  • '  f^\, ) 
is  the  vector  of  given  boundary  values.  For  the  exterior  prohlcni, 

P=if^\ fN /^f /4'm) 

and 

/  =  (/i  .  ■  •  ■ .  /n,  .  ■  ■  ■ .  /i    .  ■  ■  • .  Inm  ) 

The  N  hy  N  matrices  A''  and  A'^  in  (26)  and  (27)  represent  the  interac- 
tions of  the  discrete  double  layer  potentials.  The  A'  by  M  matrices  B'  and  B' 
represent  the  logarithmic  terms  coupling  the  density  values  to  the  coefficients 
A\,  ...,Am.  Finally,  the  A/  by  A'  matrices  C  and  C  and  the  A/  by  A/  matrices 
D'  and  D'  represent  the  discrete  constraint  equations  in  (24)  and  (25). 

4.1      Solution  of  the  Discrete  Systems 

The  matrix  equations  (26)  and  (27)  were  solved  iteralivoly,  using  the  gener- 
alized minimum  residual  method  GMRES  [22].  The  convergence  rale  of  this 
scheme  depends  on  the  eigenvalues  of  the  matrix  More  specifically,  real  and 
clustered  eigenvalues  are  advantageous,  enabling  convergence  in  a  small  number 
of  iterations. 

Observation  4.1:  Since  the  integral  operators  in  (17)  and  (23)  are  compact, 
they  can  be  approximated  to  any  finite  precision  by  an  operator  of  finite  rank. 
As  a  result,  the  eigenvalues  of  the  matrices  I  -f  A''  and  /  —  A''  can  be  shown  to 
be  cusymptotically  bounded  and  to  cluster  around  one.  In  the  continuous  case 
they  are  real,  despite  the  fact  that  the  matrices  are  nonsyminctric  (see  [13]). 

Definition  4.1:  Because  of  the  preceding  observation,  we  will  say  that  /  -I-  A'' 
and  /  -  A'  are  effectively  low  rank  perturbations  of  the  idenlity. 
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The  full  matrices  in  (26)  and  (27)  can  also  be  viewed  as  eirectively  low 
rank  perturbations  of  the  identity.  The  rank  of  the  perturbation  may  grow, 
however,  as  the  blocks  B,C  and  D  grow  in  dimension  (with  increasing  numbers 
of  components)  and  as  the  geometry  becomes  more  complicated  (affecting  the 
spectral  behavior  of  7\''  and  A''). 

In  order  to  obtain  a  matrix  whose  eigenvalues  are  better  distributed  for 
iterative  methods  such  as  GMRES,  a  preconditioner  can  be  used.  Thus,  instead 
of  solving  the  linear  systems  in  (26)  and  (27),  we  solve  the  systems 

I      fl'    \"'  /   /+  A''     B'   \  f  p 
CD'  [       C         D-       [   a 


(<^.  ^:)"(Y)- 


and 


I      B'   \       f  I  -  K'     B'   \  (  p 


CD'  \       C        D' 


a 


I       B'   V'  f  -2/ 


CD'  V      0 


(29) 


It  is  easy  to  verify  that  these  preconditioned  matrices  are  effectively  low  rank 
perturbations  of  the  identity,  where  the  rank  of  the  perturbation  is  determined 
by  the  integral  operator  terms  A'  and  A'  but  is  no  longer  affected  by  the  block 
matrices  B' ,  B' ,C .C ,  D' ,  D' . 

Remark  4.2:  In  the  remainder  of  this  section,  we  will  ignore  the  distinction 
between  the  interior  and  exterior  problems  and  omit  the  superscripts  i  and  e 
from  B,  C  D  and  A' 

At  each  iteration  now,  it  is  necessary  to  solve  a  linear  system  with  the 
preconditioning  matrix,  so  it  is  important  that  such  a  system  be  muicIi  easier 
to  solve  than  the  original  full  matrix  equation.  This  is  accomplished  as  follows. 
First  we  form  the  M  by  A/  Schur  complement  S  of  D  in  the  preconditioner 

I      B 

given  by 

S  =   D  -  CB  .  (30) 

We  then  factor  5  using  Gaussian  elimination    To  solve  the  linear  syslcm 

/    B  \  f  r,  \       f  f 


c   Dj  K,.j  -  yf.  "  (^^> 

wc  then  backsolve  with  the  triangular  factors  to  obtain  zZ  from 
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Sz'a    -    rl    -   CrJ,  (32) 

and  use  the  computed  value  of  zZ  to  obtain  zj,  from 

z-    =    r-    -   BzZ.  (33) 

The  initial  factorization  requires  ^M^  operations,  while  backsolving  requires 
M^  operations  at  each  GMRES  iteration. 

The  bulk  of  the  work  at  each  iteration,  however,  lies  in  applying  the  full 
matrix  to  a  vector.  The  product 

C         D  J    [S 

can  be  computed  in  time  0(N  +  M)  using  the  adaptive  fast  nuillipole  method  [5], 
a  recently  developed  technique  for  the  evaluation  of  potential  fields.  In  brief, 
the  dense  matrix  vector  product  A'/T  corresponds  to  evaluating  the  pairwise 
interactions  of  a  collection  of  dipoles  located  at  the  points  Qj  ,  while  the  dense 
matrix  vector  product  Ba  corresponds  to  evaluating  the  (icld  at  the  points 
Qj  due  to  a  collection  of  M  charges  located  at  the  points  Si-.  It  should  be 
mentioned  that  we  neither  form  nor  store  the  matrix  A',  so  that  the  total  storage 
requirements  of  the  algorithm  are  0{N  +  A/^),  the  latter  term  required  for  the 
preconditioning  step. 

5      The  Dirichlet-Neumann  Map  and  the  Neu- 
mann Problem 

Our  approach  to  the  Dirichlet  problem  leaves  us  with  a  representation  of  the 
solution  U(P)  in  the  form  of  a  double  layer  potential  coTubined  with  a  number 
of  logarithmic  sources  (equations  (13)  and  (18)).  The  contribution  to  -^  from 
the  logarithmic  sources  can  clearly  be  computed  explicitly,  while  ^  applied  to 
the  double  layer  potential  results  in  a  hypersingular  integral,  in  two  dimensions, 
the  theory  of  holomorphic  functions  can  be  used  to  simplify  tiiis  calculation.  We 
begin  by  observing  thai 


0{P)  =  ^  I  ti{Q)^\n\Q-P\dQ 


d_ 

'  dUQ 


is  the  real  part  of  the  Cauchy  integral 
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where  we  have  equated  the  points  P  and  Q  in  the  plane  with  the  corresponding 
complex  numbers  z  and  (.  U  and  its  harmonic  conjugate  V  satisfy  the  Cauchy- 
Riemann  equations,  so  that 

dU__dV_ 

du        dr 

and  it  remains  only  to  compute  V  and  its  tangential  derivative.  For  this  purpose, 
we  consider  the  boundary  mesh  to  be  divided  into  odd  and  even  points.  The 
value  of  V  at  the  odd  points  is  then  obtained  by  using  the  trapezoidal  rule 
applied  to  the  even  points,  and  vice-versa.  It  can  be  shown  [24]  that  such 
a  rule  is  superalgebraically  converging  on  smooth  domains.  It  has  been  used 
previously  by  a  number  of  authors  (see,  for  example,  [16], [17], [23])  Once  V  has 
been  tabulated,  we  obtain  its  tangential  derivative  via  Fourier  aj)proxiiiiation 
in  arclength  and  the  FFT. 

5.1      The  Neumann  Problem 

Many  important  problems  governed  by  Laplace's  equation  require  ijie  solution  of 
the  Neumann  problem  rather  than  the  Dirichlet  problem.  The  potential  theory 
approach  is  ecisily  applied  in  this  case  to  both  simply  connected  and  multiply 
connected  domains.  There  is  no  large  finite-dimensional  nullspace  wiiicli  requires 
special  consideration.  We  will  limit  our  attention,  for  the  sake  of  brevity,  to  the 
exterior  problem.  The  solution  of  the  Neumann  problem 


(34) 


is  sought  in  terms  of  a  single  layer  potential 

U{P)^^  [  ^iQ)^^\Q-P\dQ.  (35) 

In  order  to  satisfy  the  boundary  condition  at  Qo,  t'he  source  density  cr{Q)  must 
satisfy  the  integral  equation 


AU{P)     = 

0         for  P  €  £• 

lim       „           = 
;7?     9^ 

f(Q)         for  Q  i 

J  f(Q)dQ     = 

0  , 

^(Qo)  +  -   /  fi{Q)^^^  In  \Q  -  Qo\  dQ  =  2/(Qo)  .  (36) 

'"  Jl 


d_ 


This  equation  is  invertible  for  any  permissible  data  f(Q)  (sec  [8]) 

In   the  next   section,   we  give  one  example  of  the  solution  of  a  Neumann 
problem,  corresponding  to  potential  flow. 
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6      Numerical  Results 

We  illustrate  the  performance  of  the  algorithm  by  considering  both  interior  and 
exterior  problems  on  a  variety  of  domains. 

Example  1:  M  -  6.  The  Dirichlet  problem  is  solved  in  a  domain  containing 
six  elliptical  holes  (Fig.  1).  The  precise  centers,  eccentricities  and  inclination 
angles  of  the  ellipses  are  given  in  Table  1  and  the  data  is  obtained  by  choosing 
an  exact  solution  of  the  form 


U{x,y)^C^D 


X  -  Xo 


[r  -  xor  +  {y  -yo 


-  .mIS 


]+^/liln((j-x,)-  +  (y-j/t)2)   (37) 


where  (xo,yo)  lies  outside  the  outer  boundary  and  the  points  {x,,y,),  i  = 
1,  .  .  .,6  lie  inside  the  six  holes.  The  latter  points  are  distinct  from  the  points 
Sk  used  in  our  method.  For  this  example  we  choose 

C=  1.0,      Ak  =  /t-7/2,    k=  1 6, 

D  —  2.0  for  the  interior  problem,  and  D  =  0.0  for  the  exterior  problem. 

Table  1:  Domain  in  Example  1.  The  first  six  rows  correspond  to  the  holes,  while 
the  last  row  corresponds  to  the  outer  boundary  for  the  interior  problem. 


^(</')  -  Cr  +  acos^cosfll)  -  fcsin^sin«i>,    0  <=  (p  <=  'Iv 
y(<i>)  =  Cy  +  6cos0sin(/>  +  asin^cos«i,    0  <=  </>  <=  2;r 

a 

b 

Cr 

cy 

0 

.3626 

.1881 

.1621 

.5940 

3.3108 

.5061 

.6053 

-1.7059 

.3423 

.5778 

.6051 

.7078 

.3577 

-.9846 

4.1087 

.7928 

.3182 

1.0000 

1.2668 

2.6138 

.3923 

.4491 

-1.9306 

-1.0663 

4.4057 

.2976 

.6132 

-.8330 

-2.1650 

5,7197 

4.0000 

3.0000 

-.5000 

-.5000 

i.OOOO 

Figures  3a  and  3b  show  the  convergence  of  the  GMRES  iterative  method, 
with  and  without  preconditioning,  for  the  interior  and  exterior  problem,  re- 
spectively. In  each  case,  the  right-hand  side  of  the  equations  was  used  as  the 
initial  guess  The  Euclidean  norm  of  the  residual  in  the  preconditioned  system, 
the  pseudoresidual  IS  plotted,  in  cases  where  the  preconditioner  was  used.  One 
hundred  discretization  points  were  used  to  represent  each  of  the  six  interior 
boundaries  as  well  as  the  outer  boundary  of  the  interior  problem.  There  is  one 
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free  parameter  in  this  method,  namely  the  number  of  direction  vectors  to  be 
saved  and  orthogonalized  against,  before  restarting  the  iteration.  In  general, 
saving  more  vectors  will  result  in  fewer  iterations,  but  the  amount  of  work  per 
iteration  and  the  amount  of  storage  will  be  incre2ised.  For  these  problems,  since 
the  matrix  vector  multiply  is  rather  expensive  (typically  about  1Q0(A'  +  M) 
operations  using  the  fast  multipole  method),  it  makes  sense,  in  terms  of  total 
work,  to  save  a  number  of  vectors.  The  numbers  beside  each  curve  in  Figure 
3  are  the  number  of  direction  vectors  saved.  As  can  be  seen  from  the  figure,  if 
enough  vectors  are  saved,  then  GMRES  converges  satisfactorily  with  or  without 
the  preconditioner.  If  fewer  direction  vectors  are  saved,  however,  the  perfor- 
mance of  the  unpreconditioned  algorithm  is  quite  poor.  This  is  approximately 
the  behavior  one  would  expect  from  the  considerations  outlined  immediately 
after  Observation  4.1. 

The  eigenvalues  of  the  preconditioned  and  unpreconditioned  matrices  for  the 
interior  and  exterior  problems  are  plotted  in  Figure  4.  Note  that  in  each  case, 
most  of  the  eigenvalues  are  tightly  clustered  around  1.  Note  also,  however,  that 
the  imaginary  parts  of  the  eigenvalues  of  the  preconditioned  systems  are  much 
smaller  than  those  of  the  unpreconditioned  systems.  This  phenomenon  has  been 
observed  in  a  wide  variety  of  test  problems,  and  we  believe  that,  in  the  limit,  as 
the  discretization  becomes  finer,  the  eigenvalues  of  the  preconditioned  system 
lie  on  the  real  axis.  This  tends  to  be  an  advantage  for  the  GMRES  method. 

Table  2  summarizes  the  performance  of  the  algorithm  as  we  refine  the  spa- 
tial discretization.  The  first  column  gives  the  number  of  points  used  in  the 
discretization  of  each  boundary.  The  table  lists  the  number  of  preconditioned 
GMRES  iterations  (saving  20  direction  vectors)  required  to  reduce  the  Euclidean 
norm  of  the  pseudoresidual  below  10~^,  the  CPU  times  in  seconds  required  on 
a  Sun  Sparcslation  l'  and  the  least  squares  error  in  computing  the  Dirichlet- 
Neumann  map  at  all  boundary  discretization  points.  The  solution  was  also 
computed  at  several  test  points  inside  the  domain,  and  the  error  was  found  to 
be  less  than  10"'  (the  accuracy  to  which  the  linear  equation  is  solved)  in  each 
case. 

The  following  observations  can  be  made  from  Table  2: 

1.  The  number  of  GMRES  iterations  is  independent  of  the  number  of  dis- 
cretization points  used. 

2  The  CPU  time  grows  Imearly  with  the  number  of  discretization  points. 

3  Convergence  of  the  solution  of  the  discretized  problem  to  that  of  the  con- 
tinuou.'i  problem  is  superalgebraic  With  200  points  per  body,  the  dis- 
cretization is  accurate  to  at  least  7  decimal  places. 


M<cfercnrp  to  commercial  producls  in  this  paper  does  not  constitute  recommendation  or 
endorsmcnl  by  the  National  Institute  of  Standards  and  Technology. 
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Table  2:  Performance  of  the  algorithm  in  solving  the  exterior  and  interior  Dirich- 
let  problem  m  Example  1. 


A' 

Exterior  Problem 
#  Iterations      CPU  Time 

Error 

Interior  Problem 
#  Iterations      CPU  Time 

Error 

50 

13                    23.6 

.22E-2 

13                    30.1 

.23E-2 

100 

13                    67.1 

.39E^3 

13                    58,9 

.40E-3 

200 

13                    130.1 

.76E-7 

13                    108.0 

.lOE-6 

Example  2:  M  =  200 

We  next  consider  the  region  exterior  to  200  randomly  distributed  lioles  in 
the  plane,  pictured  in  Figure  5.  (The  frame  around  the  figure  is  not  part  of  the 
domain  of  the  problem,  but  is  used  to  identify  subregions.)  Each  hole  is  a  circle, 
with  radii  varying  from  3.45  x  10"^  to  6.88  x  10"". 

This  type  of  region  is  typical  in  certain  materials  studies,  such  as  the  simu- 
lation of  particle  coarsening  durmg  Ostwald  ripening  [16].  In  this  case,  particles 
of  one  material  are  embedded  in  a  surrounding  second-phase  medium,  and  the 
problem  is  to  follow  the  evolution  of  the  interface  boundary.  The  motion  of  this 
interface  is  determined  by  computing  the  Dirichlet-Neumann  map  with  bound- 
ary condition  equal  to  curvature.  More  precisely,  each  point  on  the  boundary 
is  assumed  to  move  in  the  normal  direction  with  speed  given  by  the  Neumann 
data. 

Figures  6a  and  6b  show  closeups  of  different  sections  of  the  domain  in  Figure 
5.  The  normal  vectors  on  the  boundaries  are  plotted,  and  their  lengths  represent 
the  relative  speeds  at  which  the  different  boundary  points  move.  For  clarity, 
the  velocity  vectors  are  plotted  only  at  every  fourth  boundary  point.  Note  the 
tendency  for  the  small  bodies  in  such  a  configuration  to  shrink  very  rapidly, 
while  the  larger  bodies  grow  and  become  non-circular. 

Our  discretization  used  100  points  per  body,  resulting  in  a  malri.x  of  or- 
der 20200.  This  was  solved  with  11  GMRES  iterations  and  required  about 
54  minutes  on  a  Sun  Sparcstation.  Had  the  GMRES  method  been  u.sed  with 
a  conventional  matrix-vector  multiply  routine,  rather  than  the  fast  multipole 
method,  the  time  required  would  have  been  about  50  hours.  If  one  used  Gaus- 
sian elimination  to  solve  this  dense  20200  by  20200  linear  system  (assuming 
enough  storage  were  available  to  hold  this  matrix!),  the  time  required  would  be 
about  3.0  months  Without  the  new  integral  equation  formulation,  200  sytems 
of  dimension  20000  would  have  to  be  solved,  and  even  the  GMRES/fast  multi- 
pole  iteration  scheme  would  require  more  than  100  hours. 

Example  3:  M  —  8 

We  computed   potential  flow  in   an  exterior  domain  with  eiglil   boundary 
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components  (Fig.  7).  The  flow  was  chosen  to  be  uniformly  to  the  right  with 
velocity  Uq  at  oo  and  to  have  zero  circulation  about  each  body.  Tliis  allows  us 
to  represent  the  velocity  potential  as 

<t>(r,?/)  =  [/qt  +  <I>2(x,y)  , 

where  4>2(i,i/)  satisfies  Laplace's  equation  with  Neumann  boundary  conditions 
equal  to  — n  •  {Uo,0)  The  streamlines  were  computed  using  a  4th  order  Runge- 
Kutta  method.  We  should  mention  that  in  order  to  compute  all  streamlines 
without  the  use  of  special  quadrature  formulas  for  the  evaluation  of  singular 
integrals,  we  overresolved  the  boundary,  placing  2000  points  on  each  component, 
resulting  in  a  system  of  dimension  16,008.  Sixteen  iterations  and  28  minutes 
were  required  to  solve  the  integral  equation  to  seven  digits  of  accuracy. 

Example  4:  A/  =  6. 

In  this  numerical  experiment,  we  used  our  algorithm  to  compute  the  ca|iaci- 
tance  matrix  for  a  system  of  cylindrical  conductors  corresponding  to  the  exterior 
problem  of  Example  1.  The  entries  of  the  capacitance  matrix  C  are  defined  by 


'  27r  Jr    du 


where  t^*-' '  is  the  solution  of  the  Dirichlet  problem  with  boundary  conditions 

f/'-"(P)     =      1  for  Pel; 

f/'-"(P)     =     0         for  P€  L,,i  5^  j  . 

Gauss'  theorem  shows  that  the  constants  C,j  are  simply  equal  to  the  const anl.s 
A,  (the  strengths  of  the  logarithmic  sources)  produced  by  the  aigoritliiii  when 
used  to  solve  Laplace's  equation  with  the  electrostatic  potential  fixed  a.s  above 
Thus,  we  produce  one  column  of  the  capacitance  matrix  with  the  solution  of 
each  Dirichlet  problem 

Using  200  points  per  boundary,  the  capacitance  matrix  is  given  by 


r  = 


/ 


+  0.835 

-0.157 

-0.467 

-0.189 

-0.012 

-0.010  \ 

-0.157 

+0,701 

-0.081 

-0.104 

-0.307 

-0.053 

-0.467 

-0.081 

+  1.070 

-0.204 

-0.054 

-0.263 

-0.189 

-0.104 

-0.204 

+0.556 

-0.025 

-0.035 

-0.012 

-0.307 

-0.054 

-0.025 

+0.717 

-0.319 

-0.010 

-0.053 

-0.263 

-0.035 

-0.319 

+  0.680  ) 

Each  Dirichlet  problem  was  solved  to  seven  digits  of  accuracy,  and  the  entire 
calculation  required  less  than  15  minutes  on  a  Sparcstation  1. 


18  REFERENCES 

7      Conclusions 

We  have  presented  a  new  integral  equation  method  for  the  solution  of  llie  Dirich- 
let  problem  in  multiply  connected  domains.  This  new  formulation  has  been 
combined  with  the  fast  multipole  method  to  produce  an  algorithm  capable  of 
solvmg  Laplace's  equation  in  domains  with  hundreds  of  distinct  boundary  com- 
ponents in  a  matter  of  minutes  on  a  Sparcstation  1.  The  asymptotic  CPU  time 
requirements  of  our  scheme  are  of  the  order  0{N  +  M)  in  the  unpreconditioned 
mode  and  of  the  order  0{N  +  A/^)  in  the  preconditioned  mode,  where  A'  is  the 
total  number  of  points  in  the  boundary  discretization  and  M  is  the  number  of 
distinct  boundary  components.  When  M  =  200  and  N  =  20,000,  about  half 
the  time  is  spent  in  solving  the  Schur  complement  system  in  the  preconditioning 
step.  As  M  grows  further,  the  preconditioner,  as  implemented,  will  dominate 
the  cost  of  the  solution  process.  To  be  able  to  handle  thousands  of  boundary 
components,  an  iterative  method  could  be  used  to  solve  the  Schur  complement 
system  rather  than  Gaussian  elimination,  but  such  a  modification  lia.s  not  been 
implemented.  We  have  included  in  our  codes  the  ability  to  solve  the  Dirichlet 
problem,  to  compute  the  Dirichlet-Neumann  map,  and  to  solve  the  Neumann 
problem  We  have  also  shown  how  our  algorithm  can  be  used  for  capacitance 
calculations.  The  method  is  currently  being  used  for  the  simulation  of  Ostwald 
ripening,  a  process  described  briefly  in  Example  2  above. 

W'hile  our  analysis  and  implementation  are  for  the  two-dimensional  case, 
there  is  a  straightforward  extension  of  our  ideas  to  three  dimension.s,  which  we 
will  report  at  a  later  date. 
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Figure  1:  A  bounded  multiply  connected  domain  D  in  the  plane.  The  outer 
boundary  of  D  is  denoted  by  Lq,  while  the  interior  boundary  curves  are  denoted 
by  Li,.  .  .,Lm. 
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Figure  2:    An   unbounded   multiply  connected   donnain  E  in  the  plane.     The 
boundaries  of  the  M  holes  are  denoted  bj'  Li , . . . ,  Lm  ■ 
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Figure  3:  The  convergence  rate  of  GMRES  with  and  without  preconditioning, 
(a)  refers  to  the  interior  problem  and  (b)  to  the  exterior  problem.  For  the 
preconditioned  system,  the  Euclidean  norm  of  the  pseudorestdual  is  plotted  (the 
residual  of  the  preconditioned  matrix).  Solid  lines  refer  to  the  preconditioned 
system,  while  dashed  lines  refer  to  the  unpreconditioned  system. 
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Figure  4:  Eigenvalues  of  the  preconditioned  and  unpreconditioned  matrices  for 
the  interior  and  exterior  problems.  Note  that  the  eigenvalues  are  tightly  clus- 
tered about  1. 
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Figure  5:  A  complicated  exterior  region,  typical  of  certain  calculations  in  ma- 
terial science,  such  as  Ostwald  ripening.  There  are  200  randomly  distributed 
holes  in  the  domain. 
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Figure  6:  Closeups  of  different  sections  of  the  domain  of  figure  5.  Plotted  are  the 
normal  vectors  on  the  boundaries,  with  lengths  proportional  to  the  Neumann 
data  generated  by  the  Dirichlet-Neumann  map.  Note  the  tendency  for  small 
bodies  to  shrink  rapidly  and  for  the  larger  ones  to  grow. 
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Figure  7:  Potential  flow  without  circulation  in  an  8-ply  connected  domain. 
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